We study the production of the lightest neutralinos in the radiative process $e^+e^- \to \tilde\chi^0_1 \tilde\chi^0_1\gamma$ in supersymmetric models with grand unification. We consider models wherein the standard model gauge group $SU(3)_c \times SU(2)_L \times U(1)_Y$ is unified into the grand unified gauge groups SU(5), or SO(10). We study this process at energies that may be accessible at a future International Linear Collider. We compare and contrast the dependence of the signal cross section on the grand unified gauge group, and different representations of the grand unified gauge group, into which the standard model gauge group is unified. We carry out a comprehensive study of the radiative production process which includes higher order QED corrections in our calculations. In addition we carry out a detailed study of the background to the signal process coming from the Standard Model radiative neutrino production $e^+e^- \to \nu \bar\nu \gamma$, as well as from the radiative production of the scalar partners of the neutrinos (sneutrinos) $e^+e^- \to \tilde\nu \tilde\nu^\ast \gamma$. The latter can be a major supersymmetric background to the radiative production of neutralinos when the sneutrinos decay invisibly. It is likely that the radiative production of the lightest neutralinos may be a viable channel to study supersymmetric partners of the Standard Model particles at the first stage of a International Linear Collider, where heavier sparticles may be too heavy to be produced in pairs.